﻿using System;

namespace LBSoft.Common.DotNetGPS
{
    public class DegreeGPS
    {
        /// <summary>
        /// 地球半径(米)
        /// </summary>
        private const double EARTH_RADIUS = 6378137.0;
        /// <summary>
        /// 角度数转换为弧度公式
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        private static double radians(double d)
        {
            return d * Math.PI / 180.0;
        }
        /// <summary>
        /// 弧度转换为角度数公式
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        private static double degrees(double d)
        {
            return d * (180 / Math.PI);
        }

        static double DEF_PI = 3.14159265359; // PI

        static double DEF_2PI = 6.28318530712; // 2*PI

        static double DEF_PI180 = 0.01745329252; // PI/180.0

        static double DEF_R = 6370693.5; // radius of earth

        /// <summary>
        /// 计算两个经纬度之间的直接距离--适用于近距离
        /// </summary>
        public static double GetDistance(Degree Degree1, Degree Degree2)
        {

            double ew1, ns1, ew2, ns2;
            double dx, dy, dew;
            double distance;

            // 角度转换为弧度

            ew1 = Degree1.X * DEF_PI180;

            ns1 = Degree1.Y * DEF_PI180;

            ew2 = Degree2.X * DEF_PI180;

            ns2 = Degree2.Y * DEF_PI180;

            // 经度差

            dew = ew1 - ew2;

            // 若跨东经和西经180 度，进行调整

            if (dew > DEF_PI)

                dew = DEF_2PI - dew;

            else if (dew < -DEF_PI)

                dew = DEF_2PI + dew;

            dx = DEF_R * Math.Cos(ns1) * dew; // 东西方向长度(在纬度圈上的投影长度)

            dy = DEF_R * (ns1 - ns2); // 南北方向长度(在经度圈上的投影长度)

            // 勾股定理求斜边长

            distance = Math.Sqrt(dx * dx + dy * dy);

            return distance;
        }



        //适用于远距离
        public static double GetLongDistance(double lon1, double lat1, double lon2, double lat2)
        {

            double ew1, ns1, ew2, ns2;

            double distance;

            // 角度转换为弧度

            ew1 = lon1 * DEF_PI180;

            ns1 = lat1 * DEF_PI180;

            ew2 = lon2 * DEF_PI180;

            ns2 = lat2 * DEF_PI180;

            // 求大圆劣弧与球心所夹的角(弧度)

            distance = Math.Sin(ns1) * Math.Sin(ns2) + Math.Cos(ns1) * Math.Cos(ns2) * Math.Cos(ew1 - ew2);

            // 调整到[-1..1]范围内，避免溢出

            if (distance > 1.0)

                distance = 1.0;

            else if (distance < -1.0)

                distance = -1.0;

            // 求大圆劣弧长度

            distance = DEF_R * Math.Acos(distance);

            return distance;

        }
        /// <summary>
        /// 计算两个经纬度之间的直接距离(google 算法)
        /// </summary>
        public static double GetDistanceGoogle(Degree Degree1, Degree Degree2)
        {
            double radLat1 = radians(Degree1.X);
            double radLng1 = radians(Degree1.Y);
            double radLat2 = radians(Degree2.X);
            double radLng2 = radians(Degree2.Y);
            double s = Math.Acos(Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Cos(radLng1 - radLng2) + Math.Sin(radLat1) * Math.Sin(radLat2));
            s = s * EARTH_RADIUS;
            s = Math.Round(s * 10000) / 10000;
            return s;
        }
        /// <summary>
        /// 以一个经纬度为中心计算出四个顶点
        /// <para>左上，左下，右上，右下</para>
        /// </summary>
        /// <param name="Degree1">坐标</param>
        /// <param name="distance">半径(米)</param>
        /// <returns></returns>
        public static Degree[] GetDegreeCoordinates(Degree Degree1, double distance)
        {
            double dlng = 2 * Math.Asin(Math.Sin(distance / (2 * EARTH_RADIUS)) / Math.Cos(Degree1.X));
            dlng = Math.Abs(degrees(dlng));//一定转换成角度数  
            double dlat = distance / EARTH_RADIUS;
            dlat = Math.Abs(degrees(dlat));//一定转换成角度数
            return new Degree[] { new Degree(Math.Round(Degree1.X - dlng,6), Math.Round(Degree1.Y + dlat,6)),//left-top
                                  new Degree(Math.Round(Degree1.X - dlng,6), Math.Round(Degree1.Y - dlat,6)),//left-bottom
                                  new Degree(Math.Round(Degree1.X + dlng,6), Math.Round(Degree1.Y + dlat,6)),//right-top
                                  new Degree(Math.Round(Degree1.X + dlng,6), Math.Round(Degree1.Y - dlat,6)) //right-bottom
            };
            //return new Degree[] { new Degree(Math.Round(Degree1.X + dlat,6), Math.Round(Degree1.Y - dlng,6)),//left-top
            //                      new Degree(Math.Round(Degree1.X - dlat,6), Math.Round(Degree1.Y - dlng,6)),//left-bottom
            //                      new Degree(Math.Round(Degree1.X + dlat,6), Math.Round(Degree1.Y + dlng,6)),//right-top
            //                      new Degree(Math.Round(Degree1.X - dlat,6), Math.Round(Degree1.Y + dlng,6)) //right-bottom
            //};
        }
    }
    /// <summary>
    /// 经纬度坐标
    /// </summary>
    public class Degree
    {
        /// <summary>
        /// 坐标
        /// </summary>
        /// <param name="x">经度</param>
        /// <param name="y">纬度</param>
        public Degree(double x, double y)
        {
            X = x;
            Y = y;
        }
        private double x;
        /// <summary>
        /// 经度：longitude
        /// </summary>
        public double X
        {
            get { return x; }
            set { x = value; }
        }
        private double y;
        /// <summary>
        /// 纬度：latitude
        /// </summary>
        public double Y
        {
            get { return y; }
            set { y = value; }
        }
    }


    public class GPSChange
    {
        private const double pi = 3.14159265358979324;
        private const double x_pi = 3.14159265358979324 * 3000.0 / 180.0;

        //克拉索天斯基椭球体参数值
        private const double a = 6378245.0;
        //第一偏心率
        private const double ee = 0.00669342162296594323;


        /// <summary>
        /// GCJ-02转换BD-09
        /// </summary>
        /// <param name="gg_lat">纬度</param>
        /// <param name="gg_lon">经度</param>
        /// <returns></returns>
        public static GPSPoint GCJ02_to_BD09(double gg_lat, double gg_lon)
        {
            GPSPoint point = new GPSPoint();
            double x = gg_lon, y = gg_lat;
            double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * x_pi);
            double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * x_pi);
            double bd_lon = z * Math.Cos(theta) + 0.0065;
            double bd_lat = z * Math.Sin(theta) + 0.006;
            point.SetLat(bd_lat);
            point.SetLng(bd_lon);
            return point;
        }

        /// <summary>
        /// BD-09转换GCJ-02
        /// </summary>
        /// <param name="bd_lat">纬度</param>
        /// <param name="bd_lon">经度</param>
        /// <returns></returns>
        public static GPSPoint BD09_to_GCJ02(double bd_lat, double bd_lon)
        {
            GPSPoint point = new GPSPoint();
            double x = bd_lon - 0.0065, y = bd_lat - 0.006;
            double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * x_pi);
            double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * x_pi);
            double gg_lon = z * Math.Cos(theta);
            double gg_lat = z * Math.Sin(theta);
            point.SetLat(gg_lat);
            point.SetLng(gg_lon);
            return point;
        }



        /// <summary>
        /// WGS-84转换GCJ-02
        /// </summary>
        /// <param name="wgLat">纬度</param>
        /// <param name="wgLon">经度</param>
        /// <returns></returns>
        public static GPSPoint WGS84_to_GCJ02(double wgLat, double wgLon)
        {
            GPSPoint point = new GPSPoint();
            if (OutOfChina(wgLat, wgLon))
            {
                point.SetLat(wgLat);
                point.SetLng(wgLon);
                return point;
            }
            double dLat = TransformLat(wgLon - 105.0, wgLat - 35.0);
            double dLon = TransformLon(wgLon - 105.0, wgLat - 35.0);
            double radLat = wgLat / 180.0 * pi;
            double magic = Math.Sin(radLat);
            magic = 1 - ee * magic * magic;
            double sqrtMagic = Math.Sqrt(magic);
            dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
            dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);
            double lat = wgLat + dLat;
            double lon = wgLon + dLon;
            point.SetLat(lat);
            point.SetLng(lon);
            return point;
        }


        public static void Transform(double wgLat, double wgLon, double[] latlng)
        {
            if (OutOfChina(wgLat, wgLon))
            {
                latlng[0] = wgLat;
                latlng[1] = wgLon;
                return;
            }
            double dLat = TransformLat(wgLon - 105.0, wgLat - 35.0);
            double dLon = TransformLon(wgLon - 105.0, wgLat - 35.0);
            double radLat = wgLat / 180.0 * pi;
            double magic = Math.Sin(radLat);
            magic = 1 - ee * magic * magic;
            double sqrtMagic = Math.Sqrt(magic);
            dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
            dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);
            latlng[0] = wgLat + dLat;
            latlng[1] = wgLon + dLon;
        }

        private static bool OutOfChina(double lat, double lon)
        {
            if (lon < 72.004 || lon > 137.8347)
                return true;
            if (lat < 0.8293 || lat > 55.8271)
                return true;
            return false;
        }

        private static double TransformLat(double x, double y)
        {
            double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0;
            ret += (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0;
            return ret;
        }

        private static double TransformLon(double x, double y)
        {
            double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(x * pi) + 40.0 * Math.Sin(x / 3.0 * pi)) * 2.0 / 3.0;
            ret += (150.0 * Math.Sin(x / 12.0 * pi) + 300.0 * Math.Sin(x / 30.0 * pi)) * 2.0 / 3.0;
            return ret;
        }

    }

    public class GPSPoint
    {
        private double lat;// 纬度
        private double lng;// 经度

        public GPSPoint()
        {
        }

        public GPSPoint(double lng, double lat)
        {
            this.lng = lng;
            this.lat = lat;
        }


        public double GetLat()
        {
            return lat;
        }
        public void SetLat(double lat)
        {
            this.lat = lat;
        }
        public double GetLng()
        {
            return lng;
        }
        public void SetLng(double lng)
        {
            this.lng = lng;
        }


    }

}
